Axion Radiation from Strings 



C. Hagmann 1 , S. Chang 2 and P. Sikivie 3 

1 Lawrence Livermore National Laboratory, Livermore, CA 94550 
2 Department of Physics, Purdue University, W. Lafayette, IN 41901 
^Department of Physics, University of Florida, Gainesville, FL 32611 

(February 1, 2008) 



Abstract 

This paper revisits the problem of the string decay contribution to the axion 
cosmological energy density. We show that this contribution is proportional 
to the average relative increase when axion strings decay of a certain quantity 
iVax which we define. We carry out numerical simulations of the evolution and 
decay of circular and non-circular string loops, of bent strings with ends held 
fixed, and of vortex-antivortex pairs in two dimensions. In the case of string 
loops and of vortex-antivortex pairs, iV ax decreases by approximately 20%. 
In the case of bent strings, remains constant or increases slightly. Our 
results imply that the string decay contribution to the axion energy density 
is of the same order of magnitude as the well-understood contribution from 
vacuum realignment. 
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I. INTRODUCTION 



The axion was proposed over two decades ago as a solution to the 'Strong CP Problem', 
i.e. to explain why QCD conserves the discrete symmetries P and CP in spite of the fact that 
the Standard Model as a whole violates those symmetries. It is the quasi-Nambu-Goldstone 
boson |lj associated with the spontaneous breakdown of a Upq(1) symmetry which Peccei 
and Quinn postulated [0. The properties of the axion depend mainly on one unknown 
parameter, the magnitude of the vacuum expectation value v a which breaks Upq(1). The 
mass of the axion and its couplings are inversely proportional to v a . The (zero-temperature) 
mass is given by: 

10 12 GeV , , 

m a ~6/xeV-iV . (1.1) 

N is a strictly positive integer that describes the color anomaly of Upq(1). The combination 
of parameters f a = % is usually called the " axion decay constant" . 

A priori v a , and therefore m a , is a free parameter of the theory. However, it is severely 
constrained by accelerator experiments and astrophysical arguments. Combined, the accel- 
erator and astrophysical constraints imply m a < 10 _2 eV ||. In addition, there is a lower 
limit on the axion mass from the requirement that axions do not overdose the Universe. 
The axion cosmological energy density receives contributions from vacuum realignment [|J, 
from string decay [pTfTTll, and from wall decay P|T^-H]. All three contributions increase 



with decreasing axion mass. The contribution whose size has been most controversial, and 
which is the topic of this paper, is that from string decay. 

Axion models contain strings because they have a spontaneously broken U(l) symmetry, 
namely Upq(1). The latter is a global symmetry and hence axion strings are global strings. 
The strings are only present when the axion is massless, but such is indeed the case in the 
very early universe, when the temperature is much larger than 1 GeV. For the purposes of 
this paper we may describe the dynamics of axions and axion strings by the action density: 



C = -d^d^-'-(<j ) ^-vl) 2 , (1.2) 

where = <pi + i<p2 is a complex scalar field. £ is invariant under the Upq(I) symmetry: 
<fi — ► e ia cf). The symmetry is spontaneously broken by the vacuum expectation value 

cf) = v a (1.3) 

where a is the axion field. A static axion string stretched along the z axis is the field 
configuration: 

* = v a fC-) e ie , (1.4) 

where (z,r,8) are cylindrical coordinates, 5 = is the string core size, and /(|) is a 
function which goes to zero when r — > 0, approaches one for r>5, and solves equations of 
motion which follow from Eq. (|1.2|). The energy per unit length of the string is 
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<Px4-\3e»\ 2 = wvlhi£) (!-5) 
p>s 2 5 



where L is an infra-red cutoff. The R.H.S. of Eq.( |1.5|) neglects the contribution to r from 
the string core. The string core contribution is of order 7rt> 2 , i.e. it is smaller than the 
contribution (1.5) of the field outside the string core by the factor ln(L/<5). In the situations 
of interest to us, hi(-|) ~ 67. Indeed, for axion strings in the early universe, the infra-red 
cutoff is provided by the presence of neighboring strings. L is then of order the distance 
between strings which, as we will see, is of order the horizon scale. The relevant time 
is the start of the QCD phase transition, and hence L ~ 1CT 7 sec. On the other hand, 
S- 1 ~ v a ~ 10 12 GeV. 

Eq. ( |L5| ) is valid for many but not all axion models. In general 

r^7r(|) 2 ln(^) (1.6) 

where K is a model-dependent integer defined as follows. Let a be the angle conjugate 
to the PQ charge. K is the factor by which the period in a of the manifold of vacuum 
expectation values of the model is reduced by the presence of exact continuous symmetry 
transformations, such as gauge transformations. In the PQWW model 0,[l| K = 2, whereas 
in the DFSZ |]§ and KSVZ @ models K = 1. 



Axion strings appear as topological defects in the early universe when Upq(1) becomes 
spontaneously broken by the vacuum expectation value (|1.3j), at a temperature of order v a . 
The phase transition where this happens is called the PQ phase transition. We assume in 
this paper that no inflation occurs after the PQ phase transition. If there is inflation after 
the PQ transition, the axion field gets homogenized over enormous distances and the axion 
strings are 'blown away'. In that case there is no string, nor wall, decay contribution to the 
axion cosmological energy density. 

At first the strings are stuck in the plasma |J, but soon the plasma becomes sufficiently 
dilute that the strings move freely and acquire relativistic speeds. The strings are present 
till the axion mass turns on at the QCD phase transition. The critical time is t\ defined by 
m a (T(t\))t\ = 1, where m a (T) is the temperature-dependent axion mass. One finds M 

h ~2- l(T 7 sec ( J a - ) . (1.7) 



10 12 GeV 



The corresponding temperature is: 



/l0 12 GeV\ 5 

T x ~ lGeV ( ; J . (1.8) 
All strings become the edges of Nd domain walls at t\, where Nd is the number of degenerate 



vacua of the axion model [17|. It is related to N by 



N 

N d =x (1.9) 
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where K is the same factor as appears in Eq. (|1.6| ). If Na > 2, the universe becomes domain 
wall dominated. The resulting cosmology is inconsistent with observation. If = 1, the 



walls bounded by string are unstable and decay into axions [[18]. We recently gave a detailed 



discussion of the wall decay contribution to the axion cosmological energy density [|14 



In this paper, our goal is to determine the present energy density p„ tr (to) i n axions that 
were radiated by axions strings between the PQ phase transition and time t\. In section 2, 
we analyse the problem theoretically. We express p„ r (to) i n terms of quantities £, x an d r 
which parametrize the main sources of uncertainty. Most of the past debate has focused on 
f which is a functional of the energy spectrum of axions emitted by strings. The remaining 
sections of the paper report on our estimates of f using numerical simulations. Section 3 
describes our simulations of the motion and decay of circular and non-circular string loops. 
Section 4 does the same for bent strings. Section 5 describes our simulations of the motion 
and annihilation of vortex-antivortex pairs in 2 (space) dimensions. The behaviour of the 
vortex-antivortex system can be predicted by analytical methods in the regime where the 
vortex and antivortex do not overlap. Thus the vortex-antivortex pair provides an interesting 
case study where theory and simulation can be confronted with each other. In section 6 we 
summarize our results and compare them to the previous simulations by two of us ||, to 
those of Battye and Shellard 0], and those of Yamaguchi, Kawasaki and Yokoyama fll||. 



II. THEORETICAL ANALYSIS 

As we shall see, the axions radiated by strings become non-relativistic soon after t\. Since 
each axion contributes m a to the energy today, our focus is on determining their number 
density n s ^ r {t). 

Axions are radiated by collapsing string loops and by oscillating wiggles on long strings. 
By definition, long strings stretch across the horizon. They move at relativistic speeds and 
intersect one another. When strings intersect, there is a high probability of reconnection, 
i.e. of rerouting of the topological flux [IT5]. Because of such 'intercommuting', long strings 



produce loops which then collapse freely. In view of this efficient decay mechanism, the 
average density of long strings is expected to be of order the minimum consistent with 
causality, namely one long string per horizon. Hence the energy density in long strings: 

Pst r (t)=e^^(|) 2 ^ln(^) , (2.1) 

where £ is a parameter of order one. 

The equations governing the number density nf Y (t) of axions radiated by axion strings 
are || 

^Pstr „ rj- dp 

~dT = - 2Hpst * ~ (2 ' 2) 

and 

,str 



dt a u(t) dt 

where u(t) is defined by: 



SHnf + rb " — (2.3) 
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dk dp s t r ^a ^ 



u(t) J k dt dk 

k is wavevector magnitude. pBt JT* a (t) is the rate at which energy density gets converted from 
strings to axions at time t, and dp ^ T ^ (t,k) is the spectrum of the axions produced. u(t) 
is therefore the average energy of axions radiated in string decay processes at time t. The 
term —2Hp str = +Hp str — 3Hp stT in Eq. ( |2.2| ) takes account of the fact that the Hubble 
expansion both stretches (+Hp BtT ) and dilutes (SHp stT ) long strings. Integrating Eqs. ( |2.1| 
- |2~3| ), setting H = ^, and neglecting terms of order one versus terms of order ln(|), one 
obtains 

nf (t) , (2.5) 



KHi Jt PQ t'iuj{t>) 

where tpQ is the time of the PQ transition. 

The central question is: what is the average energy u(t) of axions radiated in string 
decay processes at time tl Axions are radiated by wiggles on long strings and by collapsing 
string loops. Consider a process which starts at t in and ends at tg n , and which converts an 
amount of energy E from string to axions. t in and t^ n are both of order t. The number of 
axions radiated is 

N = J dk^(U n )± (2.6) 

where ^(tfi n ) is the energy spectrum of the field in the final state. The average energy of 
axions emitted is: 

"=N ■ (2 ' 7) 

It is useful to define the quantity || 

N^t) = J dk ^(t)l ■ ( 2 -8) 

The final value of iV ax is N. The initial value of A^ax is determined in terms of E by the 
fact that the energy stored in string has spectrum ^| (t- m ) ~ |. This spectral shape - equal 
energy per unit logarithmic wavevector interval - is implied by Eq. (1.5). If £ = — is the 
length of string converted to axions, we have 

f ftj = M|) 2 4 (2.9) 

for k min < k < k max where /c max is of order and k min of order ^f. t plays the role of L in 
Eq. (1.5). Hence 

AUtin) = f 1 dk± = - , t f, (2.10) 



K Jk min k 2 Mi) k 



''mm 



Combining Eqs. (|2.6| - |2.10|) yields 
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where r is the relative change in N ax (t) during the process in question: 

jVax(tfin) 



(2.12) 



k min is of order ^ where L is the loop size in the case of collapsing loops, and the wiggle 
wavelength in the case of bent strings. L is at most of order t but may be substantially 
smaller than that if the string network has a lot of small scale structure. To parametrize our 
ignorance in this matter, we define x such that the suitably averaged /c m i n = X^f- Combining 
Eqs. O and ( PH ) we find: 

~ V"^ > ( 2 - 13 ) 

X KH 

where f is the weighted average of r over the various processes that convert string to axions. 
Let us show that the set of all axions that were radiated between tpQ and t have spectrum 
~ rj for j < k < i — , irrespective of the shape of dp jT7 k a ■ Indeed scaling implies 



i^(t,k) = ^^(t)m)t, (2.i4) 

ere ct/e at 
where the unknown function f(u) is normalized such that 

f(u)du = l. (2.15) 



We will only assume that f(u) has appreciable support near u = 1. Eqs. (|2.1| , ^2|) imply 

^pW-^Pst r (t) , (2.16) 

where " ~" indicates, as before, that terms of order one are neglected versus terms of order 
ln(|). Since axions free-stream after they are emitted, and R ~ \fi at the relevant epoch, 
we have: 

with fc' = k(jj)2. The factor ~ accounts for the redshift and the volume expansion between 
t' and t. Using Eqs. ( |2.14| , |2.16| ) and (PTTj), one obtains 

/ ^WMlL) • (2-18) 



dk ' fc 2 t 2 Jk^/ttpQ kH5 

the integral in Eq. ( |2.18|) 
to the factor That is the promised result 



For 1 < k ^ — ^= — , the integral in Eq. ( |2.18| ) is only a slowly varying function of k compared 
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Since their momenta are of order tj" 1 at time t%, the axions become non-relativistic soon 
after they acquire mass. Therefore, the string decay contribution to the axion energy density 
today is 

pf (to) = m»nf fa) A 3 ~ rnX-^-{^-f (2.19) 
Ro x KHi R 

where ^ is the ratio of scale factors between t\ and today. For comparison, the contribution 
from vacuum realignment is || 

pT(h)^rnJf&Y . (2.20) 
In terms of the critical energy density p c = the vacuum realignment contribution is: 

Pc 3 l 10 12 GeV J 1 h } K } 

As usual, h parametrizes the present Hubble rate H = h - 100 scc k ^ pc - The contribution from 
wall decay is ITJI 



p^{t )^m a - f f (|i) 3 (2.22) 

where 7 is the average Lorentz factor of axions produced in the decay of walls bounded by 
string. In simulations, we found 7 ~ 7 for ^- ~ 500, but that 7 increases approximately 
linearly with ln(— ). Extrapolation of this behaviour to the parameter range of interest, 
ln(— ) ~ 60, yields 7 ~ 60. It suggests that the wall decay contribution is subdominant 
relative to the vacuum realignment contribution. 

The ratio of the string decay and vacuum realignment contributions is: 

pfM^irAg (2 23) 

^vac (+ \ •%/ ' 

Pa {to) X 

where we used Eq. ( |1.9| ). Each of the factors on the R.H.S. deserves discussion: 

Almost surely one needs = 1 to avoid the cosmological disaster of an axion 
domain wall dominated universe. It may be worth pointing out, however, that the domain 
wall problem can be avoided, in > 1 models, by introducing an interaction which slightly 



lowers one of the vacua with respect to the others |17| . The lowest vacuum takes over 
after some time and the walls disappear before they dominate the energy density. There is 
little room in parameter space for this to happen, but it is logically possible. It is discussed 



in detail in ref. 



£ In previous work @|n| , we set £ = 1 on the argument that the number density of long 
strings should be close to the minimum consistent with causality, i.e. of order one long string 
per horizon. Battye and Shellard |J set £ ~ 13 because this describes the density of strings 
in simulations of local string networks in an expanding universe. However, axion strings are 
global strings. Unlike global strings, local strings have all their energy located in the string 
core. Also, they cannot dissipate their energy by emitting Nambu-Goldstone radiation as 
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global strings do. For these reasons, it is not obvious that global strings are as dense in 
the early universe as local strings would be. In fact, M. Yamaguchi, M. Kawasaki and J. 
Yokoyama |TT] have done simulations of global string networks in an expanding universe and 
find £ ~ 1. 

X X an d £ are related since the average interstring distance controls both. On dimen- 
sional grounds, x ~ So, the effect of small scale structure in the axion string network 
partially cancels out in the RHS of Eq. ( |2.19| ). In previous work P,|TU|], we have assumed 
X — 1 but this could be off by a factor two or so. 

f This is the unknown on which most of the past debate has focused. Two basic 
scenarios have been put forth, which we call A and B. The question is: what is the spectrum 
of axions radiated by strings? The main source is closed loops of size L ~ t. Scenario 
A postulates that a bent string or closed loop oscillates many times, with period of order 
L, before it has released its excess energy and that the spectrum of radiated axions is 
concentrated near 2z. In that case one has f ~ hi(%) — 67. Scenario B postulates that 
the bent string or closed loop releases its excess energy very quickly and that the spectrum 



of radiated axions is ~ ^ with a high frequency cutoff of order and a low frequency 
cutoff of order In scenario B, the initial and final spectra 4r of the energy stored 
in the axion field are qualitatively the same and hence f ~ 1. In scenario A, the string 
decay contribution dominates over the vacuum realignment contribution by the factor ln(y), 
whereas, in scenario B, the contributions from string decay and vacuum realignment have 
the same order of magnitude. 

Computer simulations offer a way to try and estimate f. It should be kept in mind how- 
ever that present day technology limits lattice sizes to approximately 256 3 in 3 dimensions 
(3D), and 4096 2 in 2D. Hence the ratio of loop/core size that can be investigated is limited 
to hi(-|) ~ 3.5 in 3D, and 5 in 2D, whereas hi(-|) ~ 67 in the situations of physical interest. 
It is therefore important to investigate the dependence, if any, of the results of computer 
simulations on ln(j) over the small range that can be investigated. The following three 
sections report on our simulations of circular and non-circular string loops, bent strings, 
and vortex-antivortex pairs. 



III. STRING LOOP SIMULATIONS 

We simulated the motion and decay of circular loops initially at rest, and of non-circular 
loops with angular momentum. The initial configurations are set up on large (~ 10 7 points) 
Cartesian grids, and time-evolved using the finite-difference equations derived from Eq. 
( |1 . 2[ ) . FFT spectrum analysis of the kinetic and gradient energies during the collapse yields 
A ax (t), and hence r. 



A. Circular Loops 

Because of azimuthal symmetry, circular loops can be studied in p — z space. The z-axis 
is perpendicular to the plane of the loop. By mirror symmetry, the problem can be further 
reduced to one quarter-plane. The static axion field outside the string core is || 



S 



V 

z) = —to(j>, z) (3.1) 

in the infinite volume limit, where Q is the solid angle subtended by the loop. We use as 
initial configuration the outcome of a relaxation routine starting with Eq. ( |3.1| ) outside the 
core and 

Mr) = tanh(0.58^)e ie (3.2) 
o 

within the core. Here, r is the distance to the string center, and 9 is the polar angle 
about the string. The configuration inside the core is held fixed during the relaxation. 
The relaxation and subsequent dynamical evolution are done with reflective Neumann-type 
boundary conditions. A step size dt = 0.2 was used for the time evolution. The total energy 
was conserved to better than 1%. 

We call Ro the initial loop radius. Fig. 1 shows R(t) for collapsing loops. While they 
collapse, the loops reach speeds close to the speed of light. Figs. 2a through 2e show suc- 
cessive snapshots of the string core as it speeds toward the origin. It is increasingly Lorentz 
contracted. A lattice effect is observed when the Lorentz contracted core size becomes com- 
parable to the lattice spacing. This effect consists of a "scraping" of the string core on the 
underlying grid, with dissipation of the kinetic energy of the string into high frequency axion 
radiation. We chose A small enough to avoid this phenomenon. For the large (Rq = 2400) 
2D circular loop simulations, A < 0.01 is required. 

In most cases, the loops collapse without rebound. However, for the range of parameters 
80 ^ Rq/6 < 190, a bounce occurs. It is shown in Fig. 1 (solid line) for the case Ro = 2400 
and 5 = 15.8. Figs. 2f through 21 show snapshots of the axion field configuration during 
the bounce. Note that the orientation of the string switches during the bounce, i.e. a 
left-oriented loop bounces into a right-oriented one, or vice-versa. 

Spectrum analysis of the fields was performed by expanding the gradient plus kinetic 
energy 

^kin+grad = dz 2lX pdp (^0 + \v (f) ] • V0) (3.3) 

J-Lz Jo 2 2 

using 

= a mnJo{k m p) COs(k n (z + L z )) 
mn 

^^4> = ^bmnJoi^kmp) sm(k n (z + L z )) 

mn 

V p0 = c mnJo{k m p) cos(/c n (z + L z )) (3.4) 

mn 

where the k m and k n satisfy Jo(k m L p ) = and sin(2A; n L2) = respectively. The wavevector 



magnitude of mode (m, n) is k mn = yfc^ + k\. Fig. 3 shows the power spectrum at t — 
and, after the collapse, at t = 3000. At both times, it exhibits an almost flat plateau, 
consistent with dE/dk oc 1/k. The high frequency cutoff of the spectrum is increased 
however at the later time, as might be expected since the core gets Lorentz contracted 

during the collapse. The evolution of iV ax = mn was studied for various values of Rq/S. 

mn Kmn 
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As shown in Fig. 4, we observe a marked decrease of during the collapse, of order 20% 
independent of Rq/5 in the investigated range 3.6 < ln(R /5) < 5.0. In the previous circular 
loop simulations by two of us ||, in which \n(R /5) ranges from 2.5 to 3.4, it was also found 
that r ~ 0.8. 



B. Non-circular loops 

Circular loops are a special case and their behaviour may be untypical of loops in general. 
To address this concern, we performed simulations of non-circular loops as well. The initial 
condition of the loop is given by its initial position r(s, t ) and velocity v(s,t ), where 
s parametrizes distance along the loop. Only the transverse part of v(s, t ) has physical 
meaning. We determine the initial axion field using 

v 

a{x,t ) = yQ(f){f(s,t )} (3.5) 

where Q(x){f(s, to)} is the solid angle subtended by the loop as seen from x. The axion 
field a(x, to + dt) a short time dt later is similarly determined by calculating the solid angle 
subtended by the loop located at r(s, to) + dt v(s, t ). This yields the initial time derivative 
of the axion field d(x, to). 

Scenario A postulates that an axion string behaves like a Nambu-Goto (NG) string to 
lowest order in a perturbative expansion in powers of j^ry where L is the string size. The 



most general solution to the NG equations of motion is |20| : 



r(s,t) = ^-[a(a„)+b(a + )] (3.6) 

where a± = (s ± t)/Ro, s is proportional to the energy in the string between the points 
labeled s and s = 0, and a(a) and b(a) are arbitrary functions of period 2n, which satisfy 
a' (a) 2 = b\a) 2 = 1. L = 27iR is the total proper length of the string loop, i.e. f(s + L, t) = 
f(s, t). It can be shown [^TJ that a NG string loop of length L has a motion which is periodic 
in time with period L/2, and that an initially static NG string loop collapses to a doubled-up 
line after half a period. 

A considerable literature [pl]-p4| is devoted to the problem of finding non-intersecting 



NG loop solutions. The authors are mainly motivated by issues related to the cosmic string 
scenario of large scale structure formation. However, self-intersection of string loops is also 
relevant to the question whether scenario A or B is more likely to be correct. Intercommuting 
(self-intersection with reconnection) favors scenario B because it causes loop sizes to shrink, 
and hence the average energy of radiated axions to increase. As discussed in Ref. [|], if 
the probability p of intercommuting per oscillation is larger than of order K ~ 1.4%, the 

spectrum of axions radiated by the original loop, its two daughters, four grand-daughters, 
and so on, will be qualitatively the same as in scenario B, independently of assumptions on 
the spectrum of radiation from any single loop. 

A particular two-parameter set of solutions |2T| , p2| is given by: 



Rq ( . . 1 
x(s, t) = — I (1 — a) sin o\_ + —a sin 3<r_ + sin a + 
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y(s, t) = ^ f— (1 — a) cos cr_ — -a cos 3er_ — cos if> cos 

z(s, t) = ^— 2\J a{l — a) cos a_ — sin ^ cos cr + ^ . (3.7) 

with a G (0, 1) and ^ G (— 7r,7r). In a large region of (a, if)) parameter space, the loop does 
not self-intersect pif . 

We performed 21 simulations of non-circular loops using Eqs.( |3.7| ) as initial conditions, 
on a 256 3 lattice with periodic boundary conditions. Each simulation took approximately 1 
week to run. A variety of (a, if>) and A values were chosen. Figure 5 shows the collapse of 
a non-circular loop projected onto the xz plane. Here, and in all cases tried, the loop size 
shrinks to zero in one go, without oscillation or rebound. Standard Fourier techniques were 
used for the spectrum analysis, and A ax was computed as a function of time. In all but two 
of the cases tried, A ax decreases while the loop collapses. The r values depend upon a, if) 
and A, and cover the range 0.6 to 1.07. Figure 6 shows A ax (t) for Rq = 72, a = 0.7, if) = n/2 
and A = 0.0125, 0.025, 0.05 and 0.1. Some of the largest r values (1.06, 0.94, 0.90 and 0.89 
respectively) were obtained in these simulations. In them, and in all cases where ln(^-)- 
dependence was tested, r decreases with increasing ln(-^). It appears however that the 
second derivative is decreasing and that r may reach an asymptotic value for ln(^) ~ 3. 
For a = 0.7, if) = tt/2, the asymptotic value would be of order 0.9. 

Nine simulations were done for A = 0.05 and Ro = 72 (hence ln(^) = 2.78) and a variety 
of (a, if>) values. The average of r over this set is 0.77. 



IV. BENT STRING SIMULATIONS 

We have also carried out simulations of oscillating bent strings with ends held fixed. The 
oscillation amplitude decreases as the string loses energy to axion radiation. We choose the 
z axis parallel to the direction of the string when straight. L x L y L z is the size of the box. 
Initial configurations describing static, sinusoidally shaped strings were prepared using the 
ansatz 

(p(x, y, z) = tanh(^-^) expUarctanf Vo . , 2nz A ) (4.1) 

6 \ \x-x - A sm(— ))) 



where r = y(x — xq — A sin(^)) 2 + (y — yo) 2 , A is the initial amplitude, (x ,y ) = 

(lf> if) i s ^ ne equilibrium position of the string and A = L z is the string wavelength. 
We used square boxes [L x = L y ) in all cases. 

The field outside the core was thoroughly relaxed before dynamical time evolution was 
begun. Neumann boundary conditions were imposed on the four side faces (x = 0,L X , and 
y = 0,L y ) and periodic boundary conditions on the endfaces (z = 0,L Z ). The kinetic + 
gradient energy in the <f> field was spectrum analyzed at regular time intervals during the 
dynamical evolution. We took care to avoid finite volume and discrete space-time effects. 
To minimize finite volume (i.e. boundary) effects, L x ,L y > AL Z is needed. To avoid dis- 
cretization effects, A < 0.4 and time step dt < 0.2 are needed. When the latter conditions 
are satisfied, no scraping of the string on the underlying lattice is observed and total energy 
is conserved to better than one part in 10 3 . 
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We performed runs for a variety of box sizes, initial amplitudes Aq and core sizes 5 = 
Fig. 7 shows the amplitude A(t) as a function of time for initial amplitudes Aq = 20 and 
10, on a 256*256*64 lattice and A = 0.2. For t ^ 250, the damped oscillator behaviour of 
A(t) is very smooth and regular. For t <; 250, A(t) is less regular because the string is being 
driven by radiation which was emitted in the first couple of oscillations and which returned 
to the string's location after reflection by the sidewalls, 

Fig. 7 shows that, in general, A(t) is not proportional to Aq. After two oscillations the 
amplitudes are of the same order, A(t = 140) ~ 6, even though the initial amplitudes differ 
by a factor two. We confirm the following rule, already stated in ref. ||: oscillations of 
initial amplitude Aq much larger than A/10 decay rapidly, in one or two oscillations, till 
A(t) < A/10. After that the string is more weakly damped. 

Fig. 8 shows N^{t) for A = 30, 20, 10 on a 256*256*64 lattice and A = 0.2. A ax 
increases slightly, of order 1%, and then oscillates about an average value. It is not clear 
to us whether this slight increase of A ax is a real effect because we were unable to satisfy 
ourselves that it happens in the infinite volume limit. However, let us analyze it as a real 
effect. To this end, we need a formula for r in the case of bent strings, when only part of 
the string decays into axions. The fraction of string that has decayed between the initial 
time t; n and the final time t nn is given by the fractional change 1 — in the length i of 

*-in 

the string between those times. For a sinusoidally shaped string the length is determined in 
terms of the amplitude A by 

(4.2) 

Let us call N' ax the value of A ax restricted to that fraction of string which decays into axions 
between times t in and i nn . We have 

NL(tin) = (1 - ^)iV ax (t in ) (4.3) 

and 

NUten) = N ax {t Sa )-^N tx {t iu ) . (4.4) 
Hence, we estimate r for bent string simulations using the formula 

NL(Un) ^ax(tfin) - f^A ax (t in ) 




NL(tin) (1 - f-)N^(t 



(4.5) 



Table 1 shows the r values for A = 30, 20, 10 and A = 0.05, 0.1, 0.2 on a 256*256*64 
lattice, r increases with decreasing Aq, reaching values of order 1.12 for Aq = 10. However, 
for such small initial amplitudes, only a very small part of the energy stored in the string gets 
released. The larger amplitude (Aq ~ 0.5A) simulations are more relevant because they are 
typical of the cosmological setting and because proportionately more energy gets released 
in them. For large initial amplitude, r is close to one and has a tendency to decrease with 
increasing hi(-|). In contrast, scenario A predicts that r is of order ln(4). 
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V. VORTEX- ANTIVORTEX ANNIHILATION 



We also studied the annihilation of vortex-antivortex pairs in 2D. The Lagrangian is the 
same as before, Eq. (1.2). Note however that in 2D the field <fi and its expectation value v a 
have dimension of (length) -1 / 2 , and A has dimension of (length) -1 . The energy stored in a 
vortex at rest is E ~ irv \ ln(L/<5), which is analogous to Eq. (|1.5|) . The core size is 5 = 
as before. For a vortex-antivortex pair, the distance R between them is the infra-red cutoff 
L. The total energy in the pair is therefore 

E(R) ^ 2nvl\n(j) (5.1) 

when the vortex and antivortex are at rest. 

A vortex-antivortex pair has zero topological charge. It annihilates by emitting Nambu- 
Goldstone (NG) radiation. We may think of this process in 3D as the annihilation of an 
infinitely long straight string with an infinitely long parallel antistring. This is not directly 
relevant to the problem considered here since long parallel axion strings are unlikely in the 
early universe. However, because it is 2D, the process can be accurately simulated. Moreover 
the behaviour of this relatively simple system can be predicted on theoretical grounds. We 
will argue below that the energy spectrum of NG radiation emitted by a spinning vortex 
anti- vortex pair has the qualitative shape ~ 1 



dk k ' 

Simulations of vortex-antivortex annihilation were carried out previously in Refs. [|T9| , p5 



However, as far as we know, the present work is the first to spectrum analyze the NG 



radiation emitted in the decay. Ref. [26] presents simulations of the formation and evolution 
of vortices in an expanding 2+1 dimensional universe. 

What should one expect the radiation spectrum to be? It is known that, in 2+1 di- 
mensions, the axion field is dual to a gauge field which couples to the vortex as if it were 
a particle with electric charge e = \/2iiv a . This is the restriction to 2+1 dimensions of the 



well-known duality [p7p8[p5| relating the axion field in 3+1 dimensions to an anti-symmetric 



two- index gauge field A^ v {x) which couples to the world-sheet of the axion string. Hence, 
as long as the vortex and anti- vortex are at greater distance from one another than the sum 
of their core sizes (R > 25), they behave like a pair of oppositely charged particles. When 
the cores of the vortex and anti-vortex start to overlap, this description is no longer valid. 
However, in the generic situation where the pair has angular momentum, most of the energy 
has already been dissipated into radiation by then. 

The force between the vortex and anti-vortex is attractive and has magnitude 

F(R) = 2 -^ . (5.2) 

It is a manifestation of the aforementioned duality that this force can be thought of either 
as the gradient of the potential energy (571) or as the Coulomb force between two oppositely 



charged particles, of charges ±v27ru a . The acceleration of both vortex and antivortex is 
therefore: 

" """rT (5-3) 



RHf) 
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in the not highly relativistic regime, and neglecting the radiation reaction due to the emission 
of NG particles. The typical angular frequency of the radiation is u ~ a ~ ~. Hence its 
spectrum has the generic form: 

dE dEdR ,1 1 „ 1 , . 

du dRduo RuS 1 u 

i.e. it is flat on a logarithmic scale. The low frequency radiation is emitted first and the 
high frequency last. 

A 2048 2 lattice was initialized with the ansatz 

00, y) = M x -xi,y-yi) 0o(x -x 2 ,y- yi) (5.5) 

where {x\,yi) and (#2,2/2) are the locations of the vortex and antivortex respectively, and 
4>o(x,y) is the field of a single vortex at rest. A relaxation routine was applied to this 
configuration with the cores held fixed. Periodic boundary conditions were used during the 
relaxation and the subsequent dynamical evolution. The vortex and antivortex were given 
an initial relative velocity Vq perpendicular to their line of sight. The range of parameters 
simulated was < v < 0.6 for the initial velocity and 5 < 5 < 30 for the core size. Fig. 
9 shows snapshots of a decaying vortex-antivortex system for A = 0.005 and Vq = 0.5. The 
initial and final energy spectra are shown in Fig. 10. The initial spectrum, that of the 
initial vortex-antivortex pair, is flat like that of an axion string. The final spectrum is that 
of the NG radiation after the vortex and antivortex annihilated. It is also flat qualitatively, 
in agreement with the theoretical argument given above. In all cases, the final spectrum is 
found to be somewhat harder than the initial spectrum. Thus, iV ax decreases (r < 1) for all 
parameters investigated. 



VI. CONCLUSIONS 

We carried out numerical simulations of the decay of axion string into axions, for a 
variety of initial configurations. Our main goal is to estimate the factor f which appears 
in the expression for the axion cosmological energy density, Eq. ( |2.13| ). r is the relative 
increase of iV ax during a decay process, f is the average of r over the various processes that 
contribute. We simulated circular loops, non-circular loops and bent strings. 

The circular loop simulations were done in 2D, exploiting the axial symmetry. This 
allowed us to reach ln(^) ~ 5. We found that circular loops collapse in one go, except 
in the range 80 < Rq/8 < 190 where they bounce once. As far as we know, the bounce 
phenomenon was not observed in previous simulations nor was it anticipated in theoretical 
investigations. We found r ~ 0.8 for circular loops, whether or not they bounce. We did not 
find any dependence of r on ln(^-) over the range, 3.6 < ln(^ !i ) < 5.0, of the simulations. 
Our results are consistent with the previous simulations of circular loops in 3D by two of us 
||, which showed r ~ 0.8 over the range 2.6 < ln(-^) < 3.2 . 

Twenty-one non-circular loop simulations were carried out using as initial conditions a 
family of Kibble- Turok configurations parametrized by a and ip. In the case of Nambu-Goto 
strings, such initial conditions yield periodic non self-intersecting motion for most of the 
parameter space. In the simulations, the loops collapse in one go, without oscillation or 
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bounce. iV ax decreases in almost all cases, r depends on a, ip and ln(-^), and ranged from 
0.60 to 1.07. r decreases with increasing; ln(^) but appears to be reaching a limiting value 
for fn(^) ~ 3. Nine simulations were done for ln(-^) = 2.78. The average of r over this set 
is 0.77. 

The bent string simulations were done on lattices of size 256*256*64 and 256*256*128 
with the string in the direction of the shortest dimension. We found that oscillations of 
initial amplitude Aq much larger than A/10, where A is the wavelength of the wiggle, decay 
rapidly, in one or two oscillations, till A(t) ^ A/10. After that, the string is more weakly 
damped. This is consistent with what was found in ref. ||. We find that N ax increases 
slightly, by an amount of order 1%. The r values are listed in Table 1 for a representative 
set of parameter values. 

Our simulations are inconsistent with scenario A which predicts that r is of order ln(j). 
We find that r is of order one and also that r does not increase with increasing ln(-|) over 
the range investigated, 3 < hi(-|) < 5. Our bent string simulations are also inconsistent 
with scenario A in that A(t) is not proportional to the initial amplitude A . 

Battye and Shellard || have carried out bent string simulations. Their string configu- 
ration is similar to ours. They conclude that the string emits radiation mainly at twice the 
oscillation frequency of the string and that the spectrum falls off exponentially for large k. 
We do not find evidence of this. To obtain the spectrum of radiated axions they subtract 
the field of a static straight string when the string is going through the equilibrium position. 
The subtracted field is assumed to be that of the radiated axions and is spectrum analyzed. 
It is unclear to us that this subtraction procedure is valid because it neglects retardation 
and Lorentz contraction effects associated with the fact that the string is moving when it is 
going through its equilibrium position. 

Yamaguchi, Kawasaki and Yokoyama |llj] carried out simulations of a network of global 



strings in an expanding universe. They found £ ~ 1 for the parameter which describes 
the average density of strings [Eq. ([2.1Q1. We use their result below. They find that the 
spectrum of radiated axions is softer than the 4§ ~ 4 spectrum of the energy stored in 
strings, whereas we find in most cases that it is somewhat harder. It must be kept in mind 
however that, because the Yamaguchi et al. simulations describe a network of many strings, 
their effective value of hi(-|) is small compared to the range of ln(-|) values explored by the 
simulations described here. 

To estimate the axion cosmo logical energy density, we use f ~ 0.8 as implied by our 
circular and non-circular loop simulations. The contribution of bent long strings is expected 
to be much less important than that from closed loops. We use £ ~ 1 based on the simula- 



tions of global string networks in an expanding universe by Yamaguchi et al. [lT|, and set 
N d = 1 and x = 2 ±x as discussed in the Introduction. Hence our estimate for the string; 
decay contribution: 



7 



£ r ovac _ n 07 oil / fa \ 6 f0.7 



n — 7"— 0272 liraj It J • < 61) 

For the reasons mentioned in the Introduction, the wall decay contribution is probably 
subdominant compared to the vacuum realignement and string decay contributions. Hence 
our estimate for the total density in cold axions today: 
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Q a ~ f]™ c + fif = (0.4 to 0.9) 



fa 



10 12 GeV 



07 
~h 



(6.2) 



This is relevant to the dark matter searches presently in progress [30,31 , insofar that it 



suggests the range of axion masses for which axions are the dark matter of the universe. 
Let us remind the reader that the above estimate is for the case where there is no inflation 
after the PQ phase transition. It also assumes that the axion to entropy ratio is constant 
from time t± till the present. Various ways in which this assumption may be violated are 



discussed in the papers of ref. [32 



Finally, we simulated the annihilation of vortex- antivortex pairs in 2D. We find that the 
spectrum of radiated axions has the qualitative shape 4§ "~ r- This is consistent with our 
theoretical analysis of the system. We find that iV ax decreases in this case too. 
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simulations on a 256*256*64 lattice. 
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FIGURES 
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FIG. 1. Radius versus time of a collapsing circular loop, for A = 0.001 (dotted line) and 
A = 0.004 (solid line). 
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FIG. 2. Snapshots of a collapsing circular loop. The arrows give the direction of (</>i,<fe) 
internal space. The grey-scale shows the magnitude of |</>|, white for \§\ = 1 and black for \(f>\ = 



21 




FIG. 3. Spectrum of decaying circular loop for Ro = 2400, A = 0.004, at times t = (open 
boxes) and t = 5000 (closed boxes). 
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FIG. 4. iV ax versus time of a circular loop for R Q = 2400, L p = 4000, L z = 4096 and A = 0.004 
(solid), 0.001 (dashed), and 0.00025 (dotted). 
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FIG. 5. Snapshots of a collapsing non-circular loop. The simulation is carried out on a 256 3 
lattice with A = 0.05. The initial conditions are given by Eq. Q3.7| ) with a = 0.7, ip = vr/2, Rq = 72. 
The position of the string core is projeced onto the xz plane. 
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FIG. 6. 



-^ax versus time for four non-circular loop simulations on a 256 3 lattice 



conditions are given by Eq. (3.7) with Rq = 72, a 
(dot-dashed), 0.025 (dashed) and 0.0125 (solid). 
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FIG. 7. String core position of oscillating bent strings versus time. The simulation is carried 
out on a L x L y L z = 256 2 * 64 lattice with A = 64, A = 0.2,and A = 20 and 10. 
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FIG. 8. N ax for an oscillating bent string versus time. Here L x L y L z = 256 2 * 64, A = 64, 
A = 0.2 and A = 30(dotted), 20(dashed), lO(solid) 
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FIG. 9. Snapshots of collapsing vortex-antivortex pairs. The simulation is carried out on 
a 2048 2 lattice with A = 0.005 and v$ = 0.5. The initial separation between the vortex and 
anti-vortex is 1000 lattice units. 
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FIG. 10. Spectra of a vortex-antivortex pair at time times t = (open boxes) and t = 2000 
(closed boxes). The parameter values are the same as in Fig. 9. 
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